﻿//[]-----------------------------------------------------------------------------------------------------------------[]
// | 开源程序: 使用SCE-UA算法对新安江水文模型的参数进行优化 （C#语言）
// |-----------------------------------------------------------------------------------------------------------------
// | 参考文献：
// | [0]河海大学水文水资源学院 卢家波 2021 https://github.com/lujiabo98/XAJ-SCEUA
// | [1]段青云,SCEUA的原始Fortran代码,1992, https://shxy.hhu.edu.cn/2019/0904/c12296a195177/page.htm
// | [2]L. Shawn Matott改编的C++代码,2009, https://github.com/MESH-Model/MESH_Project_Baker_Creek/blob/7e0a7e588213916deb2b6c11589df0d132d9b310/Model/Ostrich/SCEUA.h
// | [3]Van Hoey S改编的Python代码,2011
// | [4]Mostapha Kalami Heris, Shuffled Complex Evolution in MATLAB (URL: https://yarpiz.com/80/ypea110-shuffled-complex-evolution), Yarpiz, 2015.
// | [5]Duan, Q.Y., Gupta, V.K. & Sorooshian, S. Shuffled complex evolution approach for effective and efficient global minimization. J Optim Theory Appl 76, 501–521 (1993). https://doi.org/10.1007/BF00939380.
// | [6]Duan, Q., Sorooshian, S., and Gupta, V. (1992), Effective and efficient global optimization for conceptual rainfall-runoff models, Water Resour. Res., 28( 4), 1015– 1031, https://doi.org/10.1029/91WR02985.
// | [7]Duan, Q., Sorooshian, S., & Gupta, V. K. (1994). Optimal use of the SCE-UA global optimization method for calibrating watershed models. Journal of Hydrology, 158(3-4), 265-284. https://doi.org/10.1016/0022-1694(94)90057-4.
// | [8]王书功. 水文模型参数估计方法及参数估计不确定性研究[M]. 河南:黄河水利出版社,2010.(https://book.douban.com/subject/5377630/)
// | [9]王书功. 水文模型参数估计方法及参数估计不确定性研究[D]. 北京:中国科学院研究生院,2006.(https://jz.docin.com/p-87849994.html)
//[]-----------------------------------------------------------------------------------------------------------------[]

using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using System.IO;

namespace SceuaApp
{
    // 模型参数类
    public class Parameter
    {
        public void ReadFromFile(string StrPath)
        {
            string filePath = Path.Combine(StrPath, "parameter.txt");
            if (File.Exists(filePath))
            {
                using (StreamReader fin = new StreamReader(filePath))
                {
                    m_KC = double.Parse(fin.ReadLine());
                    m_UM = double.Parse(fin.ReadLine());
                    m_LM = double.Parse(fin.ReadLine());
                    m_C = double.Parse(fin.ReadLine());

                    fin.ReadLine();

                    m_WM = double.Parse(fin.ReadLine());
                    m_B = double.Parse(fin.ReadLine());
                    m_IM = double.Parse(fin.ReadLine());

                    fin.ReadLine();

                    m_SM = double.Parse(fin.ReadLine());
                    m_EX = double.Parse(fin.ReadLine());
                    m_KG = double.Parse(fin.ReadLine());
                    m_KI = double.Parse(fin.ReadLine());

                    fin.ReadLine();

                    m_CS = double.Parse(fin.ReadLine());
                    m_CI = double.Parse(fin.ReadLine());
                    m_CG = double.Parse(fin.ReadLine());
                    m_CR = double.Parse(fin.ReadLine());
                    m_KE = double.Parse(fin.ReadLine());
                    m_XE = double.Parse(fin.ReadLine());
                }
            }
        }

        public void SetValues(double KC, double UM, double LM, double C,
            double WM, double B, double IM,
            double SM, double EX, double KG, double KI,
            double CS, double CI, double CG, double CR, double KE, double XE)
        {
            m_KC = KC;
            m_UM = UM;
            m_LM = LM;
            m_C = C;

            m_WM = WM;
            m_B = B;
            m_IM = IM;

            m_SM = SM;
            m_EX = EX;
            m_KG = KG;
            m_KI = KI;

            m_CS = CS;
            m_CI = CI;
            m_CG = CG;
            m_CR = CR;
            m_KE = KE;
            m_XE = XE;
        }

        public Parameter(double KC = 0.8, double UM = 20.0, double LM = 80.0, double C = 0.15,
            double WM = 120.0, double B = 0.2, double IM = 0.01,
            double SM = 20.0, double EX = 1.5, double KG = 0.35, double KI = 0.35,
            double CS = 0.1, double CI = 0.6, double CG = 0.95, double CR = 0.2, double KE = 24.0, double XE = 0.4)
        {
            SetValues(KC, UM, LM, C, WM, B, IM, SM, EX, KG, KI, CS, CI, CG, CR, KE, XE);
        }

        // 蒸散发计算相关参数
        public double m_KC;   // 流域蒸散发折算系数，敏感
        public double m_UM;   // 上层张力水容量/mm，敏感，10~50
        public double m_LM;   // 下层张力水容量/mm，敏感，60~90
        public double m_C;    // 深层蒸散发折算系数，不敏感，0.10~0.20

        // 产流计算相关参数
        public double m_WM;   // 流域平均张力水容量/mm，不敏感，120~200
        public double m_B;    // 张力水蓄水容量曲线方次，不敏感，0.1~0.4
        public double m_IM;   // 不透水面积占全流域面积的比例，不敏感

        // 水源划分相关参数
        public double m_SM;   // 表层自由水蓄水容量/mm ，敏感
        public double m_EX;   // 表层自由水蓄水容量方次，不敏感，1.0~1.5
        public double m_KG;   // 表层自由水蓄水库对地下水的日出流系数，敏感
        public double m_KI;   // 表层自由水蓄水库对壤中流的日出流系数，敏感

        // 汇流计算相关参数
        public double m_CS;   // 日模型地面径流消退系数，敏感
        public double m_CI;   // 日模型壤中流蓄水库的消退系数，敏感
        public double m_CG;   // 日模型地下水蓄水库的消退系数，敏感
        public double m_CR;   // 日模型河网蓄水消退系数，敏感
        public double m_KE;   // 马斯京根法演算参数/h，敏感，KE = N * ∆t，N为河道分段数
        public double m_XE;   // 马斯京根法演算参数，敏感，0.0~0.5
    }

    // 单个单元流域状态类
    public class State
    {
        public void SetInput(int nt, int nw, Watershed watershed)
        {
            m_P = watershed.GetP(nt, nw);
            m_EM = watershed.GetEM(nt, nw);
            m_F = watershed.GetF(nw);
        }

        public void ReadFromFile(string StrPath)
        {
            string filePath = Path.Combine(StrPath, "time.txt");
            if (File.Exists(filePath))
            {
                using (StreamReader fin = new StreamReader(filePath))
                {
                    m_dt = double.Parse(fin.ReadLine());
                }
            }
        }

        public void SetValues(double P, double EM, double F, double dt, double S0, double FR,
            double EU, double EL, double ED, double E,
            double WU, double WL, double WD, double W,
            double RIM, double R, double RS, double RI, double RG,
            double PE, double QS, double QI, double QG, double QU,
            double QU0, double[] O, double O2, double EP, double Q)
        {
            m_P = P;
            m_EM = EM;
            m_F = F;
            m_dt = dt;

            m_S0 = S0;
            m_FR = FR;
            m_EU = EU;
            m_EL = EL;
            m_ED = ED;
            m_E = E;
            m_WU = WU;
            m_WL = WL;
            m_WD = WD;
            m_W = W;
            m_RIM = RIM;
            m_R = R;
            m_RS = RS;
            m_RI = RI;
            m_RG = RG;
            m_PE = PE;
            m_QS = QS;
            m_QI = QI;
            m_QG = QG;
            m_QU = QU;
            m_QU0 = QU0;
            m_O = O;
            m_O2 = O2;
            m_EP = EP;

            m_Q = Q;
        }

        public State(double P = 0, double EM = 0, double F = 0, double dt = 24,
            double S0 = 0, double FR = 0,
            double EU = 0, double EL = 0, double ED = 0, double E = 0,
            double WU = 10, double WL = 30, double WD = 20, double W = 60,
            double RIM = 0, double R = 0, double RS = 0, double RI = 0,
            double RG = 0, double PE = 0, double QS = 0, double QI = 0,
            double QG = 0, double QU = 0, double QU0 = 0, double[] O = null,
            double O2 = 0.0, double EP = 0, double Q = 0)
        {
            SetValues(P, EM, F, dt, S0, FR, EU, EL, ED, E, WU, WL, WD, W,
                RIM, R, RS, RI, RG, PE, QS, QI, QG, QU, QU0, O, O2, EP, Q);
        }

        ~State()
        {
            m_O = null;
        }

        // 外部输入相关参数
        public double m_P;    // 单元流域降雨量，mm
        public double m_EM;   // 单元流域水面蒸发量，mm
        public double m_F;    // 单元流域面积,km2
        public double m_dt;   // 模型计算时段长,h

        // 模型状态相关参数
        public double m_S0;   // 本时段初产流面积上的平均自由水深，mm
        public double m_FR;   // 本时段产流面积比例
        public double m_EU;   // 上层蒸散发量，mm
        public double m_EL;   // 下层蒸散发量，mm
        public double m_ED;   // 深层蒸散发量，mm
        public double m_E;    // 总的蒸散发量，mm
        public double m_EP;   // 流域蒸发能力，mm
        public double m_WU;   // 上层张力水蓄量，mm
        public double m_WL;   // 下层张力水蓄量，mm
        public double m_WD;   // 深层张力水蓄量，mm
        public double m_W;    // 总的张力水蓄量，mm
        public double m_RIM;  // 不透水面积上的产流量，mm
        public double m_R;    // 总径流量，mm
        public double m_RS;   // 地面径流，mm
        public double m_RI;   // 壤中流，mm
        public double m_RG;   // 地下径流，mm
        public double m_PE;   // 净雨量，mm，PE = P - KC * EM
        public double m_QS;   // 地面径流汇流，m3/s
        public double m_QI;   // 壤中流汇流，m3/s
        public double m_QG;   // 地下径流汇流，m3/s
        public double m_QU;   // 本时段末单元流域出口流量，m3/s
        public double m_QU0;  // 上一时段末即本时段初的单元流域出口流量，m3/s
        public double[] m_O;  // 单元流域在各子河段出口断面形成的出流，m3/s
        public double m_O2;   // 单元流域在全流域出口断面形成的出流，m3/s

        // 输出外部相关参数
        public double m_Q;    // 流域出口断面流量，m3/s
    }

    // 从文本中导入降雨和蒸发数据，输出流量过程到文本中
    public class IO
    {
        public void ReadFromFile(string StrPath = "")
        {
            string pFilePath = Path.Combine(StrPath, "P.txt");
            if (File.Exists(pFilePath))
            {
                using (StreamReader fin = new StreamReader(pFilePath))
                {
                    // 读入降雨数据的行数和列数，行数表示时段数，列数表示雨量站数
                    var dimensions = fin.ReadLine().Split(' ');
                    nrows = int.Parse(dimensions[0]);
                    ncols = int.Parse(dimensions[1]);

                    // 创建二维数据，用于存储各雨量站降雨量，mm
                    m_P = new List<List<double>>(nrows);
                    for (int r = 0; r < nrows; r++)
                    {
                        m_P.Add(new List<double>(ncols));
                        var rowValues = fin.ReadLine().Split(' ');
                        for (int c = 0; c < ncols; c++)
                        {
                            m_P[r].Add(double.Parse(rowValues[c]));
                        }
                    }
                }
            }

            string emFilePath = Path.Combine(StrPath, "EM.txt");
            if (File.Exists(emFilePath))
            {
                using (StreamReader fin = new StreamReader(emFilePath))
                {
                    // 读入蒸发数据的行数和列数，行数表示时段数，列数表示蒸发站数
                    var dimensions = fin.ReadLine().Split(' ');
                    nrows = int.Parse(dimensions[0]);
                    ncols = int.Parse(dimensions[1]);

                    // 创建二维数据，用于存储各蒸发站蒸发量，mm
                    m_EM = new List<List<double>>(nrows);
                    for (int r = 0; r < nrows; r++)
                    {
                        m_EM.Add(new List<double>(ncols));
                        var rowValues = fin.ReadLine().Split(' ');
                        for (int c = 0; c < ncols; c++)
                        {
                            m_EM[r].Add(double.Parse(rowValues[c]));
                        }
                    }
                }
            }
        }

        public void WriteToFile(string StrPath = "")
        {
            string qFilePath = Path.Combine(StrPath, "Q.txt");
            using (StreamWriter fout = new StreamWriter(qFilePath))
            {
                foreach (var value in m_Q)
                {
                    fout.WriteLine(value);
                }
            }
        }

        public IO()
        {
            m_Q = null;
            nrows = 0;
            ncols = 0;
            m_P = new List<List<double>>();
            m_EM = new List<List<double>>();
        }

        ~IO()
        {
            m_Q = null;
        }

        public List<List<double>> m_P;  // 各雨量站逐时段降雨量，mm
        public List<List<double>> m_EM; // 各蒸发站逐时段水面蒸发量，mm
        public double[] m_Q;           // 流域出口断面流量过程，m3/s
        public int nrows;              // 数据行数
        public int ncols;              // 数据列数
    }

    // 流域分块
    public class Watershed
    {
        public void ReadFromFile(string StrPath = "")
        {
            string filePath = Path.Combine(StrPath, "watershed.txt");
            if (File.Exists(filePath))
            {
                using (StreamReader fin = new StreamReader(filePath))
                {
                    m_name = fin.ReadLine();
                    m_area = double.Parse(fin.ReadLine());
                    m_numRainfallStation = int.Parse(fin.ReadLine());
                    m_numEvaporationStation = int.Parse(fin.ReadLine());
                    m_numSubWatershed = int.Parse(fin.ReadLine());

                    // 读取单元流域面积
                    m_areaSubWatershed = new double[m_numSubWatershed];
                                       
                    var rowValues1 = fin.ReadLine().Split('\t');

                    for (int i = 0; i < m_numSubWatershed; i++)
                    {
                        m_areaSubWatershed[i] = double.Parse(rowValues1[i]);
                    }

                    // 读取雨量站分配比例
                    m_rateRainfallStation = new List<List<double>>(m_numSubWatershed);
                    for (int i = 0; i < m_numSubWatershed; i++)
                    {
                        m_rateRainfallStation.Add(new List<double>(m_numRainfallStation));
                        var rowValues = fin.ReadLine().Split('\t');
                        for (int j = 0; j < m_numRainfallStation; j++)
                        {
                            m_rateRainfallStation[i].Add(double.Parse(rowValues[j]));
                        }
                    }

                    // 读取蒸发站分配比例
                    m_rateEvaporationStation = new List<List<double>>(m_numSubWatershed);
                    for (int i = 0; i < m_numSubWatershed; i++)
                    {
                        m_rateEvaporationStation.Add(new List<double>(m_numEvaporationStation));
                        var rowValues = fin.ReadLine().Split('\t');
                        for (int j = 0; j < m_numEvaporationStation; j++)
                        {
                            m_rateEvaporationStation[i].Add(double.Parse(rowValues[j]));
                        }
                    }

                    // 读取雨量站点名称
                    m_nameRainfallStation = new string[m_numRainfallStation];
                    for (int j = 0; j < m_numRainfallStation; j++)
                    {
                        m_nameRainfallStation[j] = fin.ReadLine();
                    }

                    // 读取蒸发站点名称
                    m_nameEvaporationStation = new string[m_numEvaporationStation];
                    for (int j = 0; j < m_numEvaporationStation; j++)
                    {
                        m_nameEvaporationStation[j] = fin.ReadLine();
                    }
                }
            }
        }

        public void SetValues(string name, double area,
            int numRainfallStation, int numEvaporationStation, int numSubWatershed,
            double[] areaSubWatershed, List<List<double>> rateRainfallStation,
            List<List<double>> rateEvaporationStation,
            string[] nameRainfallStation, string[] nameEvaporationStation,
            List<List<double>> P, List<List<double>> EM)
        {
            m_name = name;
            m_area = area;
            m_numRainfallStation = numRainfallStation;
            m_numEvaporationStation = numEvaporationStation;
            m_numSubWatershed = numSubWatershed;
            m_areaSubWatershed = areaSubWatershed;
            m_rateRainfallStation = rateRainfallStation;
            m_rateEvaporationStation = rateEvaporationStation;
            m_nameRainfallStation = nameRainfallStation;
            m_nameEvaporationStation = nameEvaporationStation;
            m_P = P;
            m_EM = EM;
        }

        public void Calculate(IO io)
        {
            int nrows = io.nrows;   // 记录条数，即时段数
            int ncols = m_numSubWatershed;    // 单元流域个数

            // 计算各单元流域逐时段降雨量，mm
            m_P = new List<List<double>>(nrows);
            for (int r = 0; r < nrows; r++)
            {
                m_P.Add(new List<double>(ncols));
                for (int c = 0; c < ncols; c++)
                {
                    m_P[r].Add(0.0);
                    for (int i = 0; i < m_numRainfallStation; i++)
                    {
                        m_P[r][c] += io.m_P[r][i] * m_rateRainfallStation[c][i];   // 按比例计算单元流域降雨量
                    }
                }
            }

            // 计算各单元流域逐时段水面蒸发量，mm
            m_EM = new List<List<double>>(nrows);
            for (int r = 0; r < nrows; r++)
            {
                m_EM.Add(new List<double>(ncols));
                for (int c = 0; c < ncols; c++)
                {
                    m_EM[r].Add(0.0);
                    for (int i = 0; i < m_numEvaporationStation; i++)
                    {
                        m_EM[r][c] += io.m_EM[r][i] * m_rateEvaporationStation[c][i];   // 按比例计算单元流域水面蒸发量
                    }
                }
            }
        }

        public double GetP(int nt, int nw)
        {
            return m_P[nt][nw];
        }

        public double GetEM(int nt, int nw)
        {
            return m_EM[nt][nw];
        }

        public double GetF(int nw)
        {
            return m_areaSubWatershed[nw];
        }

        public int GetnW()
        {
            return m_numSubWatershed;
        }

        public Watershed()
        {
            m_name = "默认流域";
            m_area = 0.0;
            m_numRainfallStation = 0;
            m_numEvaporationStation = 0;
            m_numSubWatershed = 0;
            m_areaSubWatershed = null;
            m_rateRainfallStation = new List<List<double>>();
            m_rateEvaporationStation = new List<List<double>>();
            m_nameRainfallStation = null;
            m_nameEvaporationStation = null;
            m_P = new List<List<double>>();
            m_EM = new List<List<double>>();
        }

        ~Watershed()
        {
            m_areaSubWatershed = null;
        }

        public string m_name;                  // 流域名称
        public double m_area;                 // 流域控制面积，km2
        public int m_numRainfallStation;      // 雨量站个数
        public int m_numEvaporationStation;   // 蒸发站个数
        public int m_numSubWatershed;         // 单元流域个数
        public double[] m_areaSubWatershed;   // 单元流域面积，km2
        public List<List<double>> m_rateRainfallStation;  // 各单元流域对应各雨量站比例，按泰森多边形分配
        public List<List<double>> m_rateEvaporationStation;  // 各单元流域对应各蒸发站比例，按泰森多边形分配
        public string[] m_nameRainfallStation;  // 雨量站点名称
        public string[] m_nameEvaporationStation;  // 蒸发站点名称
        public List<List<double>> m_P;          // 各单元流域逐时段降雨量，mm
        public List<List<double>> m_EM;         // 各单元流域逐时段水面蒸发量，mm
    }

    public class Confluence
    {
        public void SetParameter(Parameter parameter) // 设置参数
        {
            CS = parameter.m_CS;
            CI = parameter.m_CI;
            CG = parameter.m_CG;
            CR = parameter.m_CR;
            IM = parameter.m_IM;
        }

        public void SetState(State state) // 设置状态
        {
            RIM = state.m_RIM;
            RS = state.m_RS;
            RI = state.m_RI;
            RG = state.m_RG;
            QS = state.m_QS;
            QI = state.m_QI;
            QG = state.m_QG;
            QU0 = state.m_QU;
            F = state.m_F;
            dt = state.m_dt;
        }

        public void UpdateState(State state) // 更新状态
        {
            state.m_QS = QS;
            state.m_QI = QI;
            state.m_QG = QG;
            state.m_QU = QU;
            state.m_QU0 = QU0;
        }

        public void Calculate() // 坡地汇流和河网汇流
        {
            // 把透水面积上的产流量均摊到单元流域上
            RS *= (1 - IM);
            RI *= (1 - IM);
            RG *= (1 - IM);

            // 出流系数换算
            M = 24.0 / dt;   // 一天划分的计算时段数
            CSD = Math.Pow(CS, 1.0 / M);   // 计算时段内地面径流蓄水库的消退系数
            CID = Math.Pow(CI, 1.0 / M);   // 计算时段内壤中流蓄水库的消退系数
            CGD = Math.Pow(CG, 1.0 / M);   // 计算时段内地下水蓄水库的消退系数
            CRD = Math.Pow(CR, 1.0 / M);   // 计算时段内河网蓄水消退系数

            // 坡地汇流
            U = F / 3.6 / 24;   // 单位转换系数
            QS = CSD * QS + (1 - CSD) * (RS + RIM) * U;   // 地面径流流入地面径流水库，经过消退(CSD)，成为地面径流对河网的总入流QS
            QI = CID * QI + (1 - CID) * RI * U;           // 壤中流流入壤中流水库，经过消退(CID)，成为壤中流对河网的总入流QI
            QG = CGD * QG + (1 - CGD) * RG * U;           // 地下径流进入地下水蓄水库，经过地下水蓄水库的消退(CGD)，成为地下水对河网的总入流QG

            QT = QS + QI + QG;

            // 河网汇流，采用线性水库法，且仅当单元流域面积大于200km2时才计算河网汇流
            if (F < 200)
            {
                QU = QT;   // 单元流域面积不大且河道较短，对水流运动的调蓄作用通常较小，将这种调蓄作用合并在地面径流和地下径流中一起考虑所带来的误差通常可以忽略
            }
            else
            {
                QU = CRD * QU + (1 - CRD) * QT;   // 线性水库法，只有在单元流域面积较大或流域坡面汇流极其复杂时才考虑单元面积内的河网汇流
            }
        }

        public Confluence(double cs = 0.1, double ci = 0.6, double cg = 0.95, double cr = 0.2, double im = 0.01,
                          double qs = 0.0, double qi = 0.0, double qg = 0.0, double qt = 0.0, double qu = 0.0,
                          double rs = 0.0, double ri = 0.0, double rg = 0.0, double rim = 0.0,
                          double qi0 = 0.0, double qg0 = 0.0, double qu0 = 0.0, double f = 0.0,
                          double u = 0.0, double m = 0.0, double csd = 0.0, double cid = 0.0, double cgd = 0.0, double crd = 0.0, double dt_ = 24)
        {
            CS = cs;
            CI = ci;
            CG = cg;
            CR = cr;
            IM = im;

            QS = qs;
            QI = qi;
            QG = qg;
            QT = qt;
            QU = qu;

            RS = rs;
            RI = ri;
            RG = rg;
            RIM = rim;

            QI0 = qi0;
            QG0 = qg0;
            QU0 = qu0;

            F = f;
            U = u;
            M = m;
            CSD = csd;
            CID = cid;
            CGD = cgd;
            CRD = crd;
            dt = dt_;
        }

        ~Confluence()
        {
        }

        // ========模型参数======== //
        public double CS;   // 地面径流消退系数，敏感
        public double CI;   // 壤中流消退系数，敏感
        public double CG;   // 地下水消退系数，敏感
        public double CR;   // 河网蓄水消退系数，敏感
        public double IM;   // 不透水面积占全流域面积的比例，不敏感

        // ========模型状态======== //
        public double QS;   // 单元流域地面径流，m3/s
        public double QI;   // 单元流域壤中流，m3/s
        public double QG;   // 单元流域地下径流，m3/s
        public double QT;   // 单元流域河网总入流（进入单元面积的地面径流、壤中流和地下径流之和），m3/s
        public double QU;   // 单元流域出口流量，m3/s

        public double RS;   // 地面径流量，mm
        public double RI;   // 壤中流径流量，mm
        public double RG;   // 地下径流量，mm
        public double RIM;  // 不透水面积上的产流量，mm

        public double QI0;  // QI(t-1)，前一时刻壤中流，m3/s
        public double QG0;  // QG(t-1)，前一时刻地下径流，m3/s
        public double QU0;  // QU(t-1)，前一时刻单元流域出口流量，m3/s

        public double F;    // 单元流域面积，km2
        public double U;    // 单位转换系数
        public double M;    // 一天划分的计算时段数
        public double CSD;  // 计算时段内地面径流蓄水库的消退系数
        public double CID;  // 计算时段内壤中流蓄水库的消退系数
        public double CGD;  // 计算时段内地下水蓄水库的消退系数
        public double CRD;  // 计算时段内河网蓄水消退系数
        public double dt;   // 模型计算时段长,h
    }

    public class Evapotranspiration
    {
        public void SetParameter(Parameter parameter) // 设置参数
        {
            KC = parameter.m_KC;
            LM = parameter.m_LM;
            C = parameter.m_C;
        }

        public void SetState(State state) // 设置状态
        {
            WU = state.m_WU;
            WL = state.m_WL;
            P = state.m_P;
            EM = state.m_EM;
        }

        public void UpdateState(State state) // 更新状态
        {
            state.m_EP = EP;
            state.m_E = E;
            state.m_EU = EU;
            state.m_EL = EL;
            state.m_ED = ED;
        }

        public void Calculate() // 计算三层蒸散发量
        {
            // 三层蒸散发计算
            EP = KC * EM; // 计算流域蒸发能力

            if (P + WU >= EP)
            {
                EU = EP;
                EL = 0;
                ED = 0;
            }
            else
            {
                EU = P + WU;

                if (WL >= C * LM)
                {
                    EL = (EP - EU) * WL / LM;
                    ED = 0;
                }
                else
                {
                    if (WL >= C * (EP - EU))
                    {
                        EL = C * (EP - EU);
                        ED = 0;
                    }
                    else
                    {
                        EL = WL;
                        ED = C * (EP - EU) - EL;
                    }
                }
            }

            // 计算总的蒸散发量
            E = EU + EL + ED;
        }

        public Evapotranspiration(double kc = 0.8, double lm = 80.0, double c = 0.15,
                                  double wu = 10, double wl = 30, double ep = 0.0,
                                  double e = 0.0, double eu = 0.0, double el = 0.0, double ed = 0.0,
                                  double p = 0.0, double em = 0.0) // 构造函数
        {
            KC = kc;
            LM = lm;
            C = c;

            WU = wu;
            WL = wl;
            EP = ep;
            E = e;
            EU = eu;
            EL = el;
            ED = ed;

            P = p;
            EM = em;
        }

        ~Evapotranspiration() // 析构函数
        {
        }

        // ========模型参数======== //
        public double KC; // 流域蒸散发折算系数，敏感
        public double LM; // 下层张力水容量/mm，敏感，60~90
        public double C;  // 深层蒸散发折算系数，不敏感，0.10~0.20

        // ========模型状态======== //
        public double WU; // 上层张力水蓄量，mm
        public double WL; // 下层张力水蓄量，mm
        public double EP; // 单元流域蒸发能力，mm
        public double E;  // 总的蒸散发量，mm
        public double EU; // 上层蒸散发量，mm
        public double EL; // 下层蒸散发量，mm
        public double ED; // 深层蒸散发量，mm

        // ========外部输入======== //
        public double P;  // 降雨量，mm
        public double EM; // 水面蒸发量，mm
    }

    public class Muskingum
    {
        public void SetParameter(Parameter parameter) // 设置参数
        {
            KE = parameter.m_KE;
            XE = parameter.m_XE;
        }

        public void SetState(State state) // 设置状态
        {
            I1 = state.m_QU0;
            I2 = state.m_QU;
            O = state.m_O;
            dt = state.m_dt;
        }

        public void UpdateState(State state) // 更新状态
        {
            state.m_O = O;
            state.m_O2 = O2;
        }

        public void Calculate() // 计算河道汇流
        {
            KL = dt; // 为了保证马斯京根法的两个线性条件，每个单元河取 KL = ∆t

            N = (int)(KE / KL); // 单元河段数

            XL = 0.5 - N * (1 - 2 * XE) / 2; // 计算单元河段XL

            double denominator = 0.5 * dt + KL - KL * XL; // 计算分母

            C0 = (0.5 * dt - KL * XL) / denominator;
            C1 = (0.5 * dt + KL * XL) / denominator;
            C2 = (-0.5 * dt + KL - KL * XL) / denominator;

            if (O == null)
            {
                O = new double[N]; // 创建存储单元流域在子河段出口断面的出流量的数组

                for (int n = 0; n < N; n++)
                {
                    O[n] = 0.0; // 单元流域在子河段出口断面的出流量为0
                }
            }

            for (int n = 0; n < N; n++)
            {
                O1 = O[n]; // 子河段时段初出流量

                O2 = C0 * I2 + C1 * I1 + C2 * O1; // 计算时段末单元流域在子河段出口断面的出流量，m3/s

                O[n] = O2; // 更新子河段时段初出流量

                I1 = O1; // 上一河段时段初出流为下一河段时段初入流

                I2 = O2; // 上一河段时段末出流为下一河段时段末入流
            }
        }

        public Muskingum(double ke = 0.0, double xe = 0.0,
                         double kl = 0.0, double xl = 0.0,
                         int n = 1, double c0 = 0.0,
                         double c1 = 0.0, double c2 = 0.0,
                         double i1 = 0.0, double i2 = 0.0,
                         double o1 = 0.0, double o2 = 0.0,
                         double[] o = null, double dt_ = 24) // 构造函数
        {
            KE = ke;
            XE = xe;
            KL = kl;
            XL = xl;
            N = n;
            C0 = c0;
            C1 = c1;
            C2 = c2;
            I1 = i1;
            I2 = i2;
            O1 = o1;
            O2 = o2;
            O = o;
            dt = dt_;
        }

        ~Muskingum() // 析构函数
        {
            // 如果需要释放资源，可以在这里处理
        }

        // ========模型参数======== //
        public double KE; // 马斯京根法演算参数，h，敏感，KE = N * ∆t
        public double XE; // 马斯京根法演算参数，敏感，0.0~0.5

        // ========模型状态======== //
        public double KL; // 子河段的马斯京根法演算参数/h
        public double XL; // 子河段的马斯京根法演算参数
        public int N;     // 单元河段数，即分段数
        public double C0; // 马斯京根流量演算公式I2系数
        public double C1; // 马斯京根流量演算公式I1系数
        public double C2; // 马斯京根流量演算公式O1系数
        public double I1; // 时段初的河段入流量
        public double I2; // 时段末的河段入流量
        public double O1; // 时段初的河段出流量
        public double O2; // 时段末的河段出流量
        public double[] O; // 各子河段出流量
        public double dt; // 模型计算时段长,h
    }

    public class Runoff
    {
        public void SetParameter(Parameter parameter) // 设置参数
        {
            WM = parameter.m_WM;
            B = parameter.m_B;
            IM = parameter.m_IM;
            WUM = parameter.m_UM;
            WLM = parameter.m_LM;
        }

        public void SetState(State state) // 设置状态
        {
            WU = state.m_WU;
            WL = state.m_WL;
            WD = state.m_WD;
            W = state.m_W;
            P = state.m_P;
            EU = state.m_EU;
            EL = state.m_EL;
            ED = state.m_ED;
            EP = state.m_EP;
        }

        public void UpdateState(State state) // 更新状态
        {
            state.m_WU = WU;
            state.m_WL = WL;
            state.m_WD = WD;
            state.m_W = W;
            state.m_R = R;
            state.m_PE = PE;
            state.m_RIM = RIM;
        }

        public void Calculate() // 产流计算并更新土壤含水量
        {
            // ========计算产流======== //
            WMM = (1 + B) / (1 - IM) * WM; // 包气带蓄水容量最大值，mm

            A = WMM * (1 - Math.Pow(1 - W / WM, 1 / (1 + B))); // 初始土壤含水量最大值，mm

            PE = P - EP;

            if (PE <= 1e-5) // 这里认为净雨量小于1e-5时即为小于等于0
            {
                R = 0.0;
                RIM = 0.0; // 计算不透水面积上的产流量
            }
            else
            {
                if (A + PE <= WMM)
                {
                    R = PE + W - WM + WM * Math.Pow(1 - (A + PE) / WMM, B + 1);
                }
                else
                {
                    R = PE - (WM - W);
                }

                RIM = PE * IM; // 计算不透水面积上的产流量
            }

            // ========计算下一时段初土壤含水量======== //
            WU = WU + P - EU - R;
            WL = WL - EL;
            WD = WD - ED;

            if (WD < 0)
            {
                WD = 0; // 防止深层张力水蓄量小于0
            }

            // 防止张力水蓄量超上限
            if (WU > WUM)
            {
                WL = WL + WU - WUM;
                WU = WUM;
            }

            if (WL > WLM)
            {
                WD = WD + WL - WLM;
                WL = WLM;
            }

            WDM = WM - WUM - WLM; // 计算深层张力水容量

            if (WD > WDM)
            {
                WD = WDM;
            }

            // 计算土壤含水量
            W = WU + WL + WD;
        }

        public Runoff(double wm = 120.0, double b = 0.2, double im = 0.01,
                      double wum = 20.0, double wlm = 80.0, double wdm = 20.0,
                      double r = 0.0, double rim = 0.0, double w = 0.0,
                      double wu = 10.0, double wl = 30.0, double wd = 20.0,
                      double wmm = 0.0, double a = 0.0,
                      double eu = 0.0, double el = 0.0, double ed = 0.0,
                      double ep = 0.0, double pe = 0.0, double p = 0.0) // 构造函数
        {
            WM = wm;
            B = b;
            IM = im;
            WUM = wum;
            WLM = wlm;
            WDM = wdm;

            R = r;
            RIM = rim;
            W = w;
            WU = wu;
            WL = wl;
            WD = wd;
            WMM = wmm;
            A = a;
            EU = eu;
            EL = el;
            ED = ed;
            EP = ep;
            PE = pe;

            P = p;
        }

        ~Runoff() // 析构函数
        {
        }

        // ========模型参数======== //
        public double WM; // 流域平均张力水容量/mm，不敏感，120~200
        public double B;  // 张力水蓄水容量曲线方次，不敏感，0.1~0.4
        public double IM; // 不透水面积占全流域面积的比例，不敏感
        public double WUM; // 上层张力水容量/mm，敏感，10~50
        public double WLM; // 下层张力水容量/mm，敏感，60~90
        public double WDM; // 深层张力水容量/mm，等于WM - WUM - WLM，不属于参数

        // ========模型状态======== //
        public double R;  // 总径流量，mm
        public double RIM; // 不透水面积上的产流量，mm
        public double W;   // 流域平均初始土壤含水量，mm
        public double WU;  // 上层张力水蓄量，mm
        public double WL;  // 下层张力水蓄量，mm
        public double WD;  // 深层张力水蓄量，mm
        public double WMM; // 包气带蓄水容量最大值，mm
        public double A;   // 初始土壤含水量最大值，mm
        public double EU;  // 上层蒸散发量，mm
        public double EL;  // 下层蒸散发量，mm
        public double ED;  // 深层蒸散发量，mm
        public double EP;  // 流域蒸发能力，mm
        public double PE;  // 净雨量，mm

        // ========外部输入======== //
        public double P;   // 降雨量，mm
    }

    public class Source
    {
        public void SetParameter(Parameter parameter) // 设置参数
        {
            SM = parameter.m_SM;
            EX = parameter.m_EX;
            KG = parameter.m_KG;
            KI = parameter.m_KI;
        }

        public void SetState(State state) // 设置状态
        {
            R = state.m_R;
            PE = state.m_PE;
            FR = state.m_FR;
            S0 = state.m_S0;
            dt = state.m_dt;
        }

        public void UpdateState(State state) // 更新状态
        {
            state.m_RS = RS;
            state.m_RI = RI;
            state.m_RG = RG;
            state.m_S0 = S;
            state.m_FR = FR;
        }

        public void Calculate() // 划分三水源
        {
            // 出流系数换算
            M = 24.0 / dt; // 一天划分的计算时段数

            KID = (1 - Math.Pow(1 - (KI + KG), 1.0 / M)) / (1 + KG / KI); // 表层自由水蓄水库对壤中流的计算时段出流系数，敏感
            KGD = KID * KG / KI; // 表层自由水蓄水库对地下水的计算时段出流系数，敏感

            // 三分水源
            if (PE <= 1e-5)
            {
                // 净雨量小于等于0时
                RS = 0;

                RI = KID * S0 * FR; /* 当净雨量小于等于0时，消耗自由水蓄水库中的水
                                  产流面积比例仍为上一时段的比例不变 */
                RG = KGD * S0 * FR;

                S = S0 * (1 - KID - KGD); // 更新下一时段初的自由水蓄量
            }
            else
            {
                // 净雨量大于0时
                SMM = SM * (1 + EX); // 全流域单点最大的自由水蓄水容量，mm

                FR0 = FR; // 上一时段产流面积比例
                FR = R / PE; // 计算本时段产流面积比例

                if (FR > 1) // 如果FR由于小数误差而计算出大于1的情况，则强制置为1
                {
                    FR = 1;
                }

                S = S0 * FR0 / FR;

                N = (int)(PE / 5.0) + 1; // N 为计算时段分段数，每一段为计算步长

                Q = PE / N; // Q 是每个计算步长内的净雨量，mm

                KIDD = (1 - Math.Pow(1 - (KID + KGD), 1.0 / N)) / (1 + KGD / KID); // 表层自由水蓄水库对壤中流的计算步长出流系数，敏感
                KGDD = KIDD * KGD / KID; // 表层自由水蓄水库对地下水的计算步长出流系数，敏感

                // 把该时段的RS、RI、RG置0，用于后续累加计算步长内的RSD、RID、RGD
                RS = 0.0;
                RI = 0.0;
                RG = 0.0;

                // 计算产流面积上最大一点的自由水蓄水容量 SMMF
                if (EX == 0.0)
                {
                    SMMF = SMM; // EX等于0时，流域自由水蓄水容量分布均匀
                }
                else
                {
                    // 假定SMMF与产流面积FR及全流域上最大点的自由水蓄水容量SMM仍为抛物线分布
                    SMMF = (1 - Math.Pow(1 - FR, 1.0 / EX)) * SMM;
                }

                SMF = SMMF / (1.0 + EX);

                // 将每个计算时段的入流量R分成N段，计算各个计算步长内的RSD、RID、RGD，再累加得到RS、RI、RG
                for (int i = 1; i <= N; i++)
                {
                    if (S > SMF)
                    {
                        S = SMF;
                    }

                    AU = SMMF * (1 - Math.Pow(1 - S / SMF, 1.0 / (1 + EX)));

                    if (Q + AU <= 0)
                    {
                        RSD = 0;
                        RID = 0;
                        RGD = 0;
                        S = 0;
                    }
                    else
                    {
                        if (Q + AU >= SMMF)
                        {
                            // 计算步长内的净雨量进入自由水蓄水库后，使得自由水蓄水深超过产流面积上最大单点自由水蓄水深
                            RSD = (Q + S - SMF) * FR;
                            RID = SMF * KIDD * FR;
                            RGD = SMF * KGDD * FR;
                            S = SMF * (1 - KIDD - KGDD);
                        }
                        else
                        {
                            // 自由水蓄水深未超过产流面积上最大单点自由水蓄水深
                            RSD = (S + Q - SMF + SMF * Math.Pow(1 - (Q + AU) / SMMF, 1 + EX)) * FR;
                            RID = (S * FR + Q * FR - RSD) * KIDD;
                            RGD = (S * FR + Q * FR - RSD) * KGDD;
                            S = S + Q - (RSD + RID + RGD) / FR;
                        }
                    }

                    // 累加计算时段内的地面径流、壤中流和地下径流
                    RS += RSD;
                    RI += RID;
                    RG += RGD;
                }
            }
        }

        public Source(double sm = 20.0, double ex = 1.5, double kg = 0.35, double ki = 0.35,
                      double r = 0.0, double rs = 0.0, double ri = 0.0, double rg = 0.0,
                      double pe = 0.0, double fr = 0.0, double s0 = 0.0, double s = 0.0,
                      double m = 0.0, double kid = 0.0, double kgd = 0.0,
                      double smm = 0.0, double smmf = 0.0, double smf = 0.0, double au = 0.0,
                      double rsd = 0.0, double rid = 0.0, double rgd = 0.0, double fr0 = 0.0,
                      int n = 1, double q = 0.0, double kidd = 0.0, double kgdd = 0.0, double dt_ = 0.0) // 构造函数
        {
            SM = sm;
            EX = ex;
            KG = kg;
            KI = ki;

            R = r;
            RS = rs;
            RI = ri;
            RG = rg;
            PE = pe;
            FR = fr;
            S0 = s0;
            S = s;

            M = m;
            KID = kid;
            KGD = kgd;

            SMM = smm;
            SMMF = smmf;
            SMF = smf;
            AU = au;

            RSD = rsd;
            RID = rid;
            RGD = rgd;
            FR0 = fr0;

            N = n;
            Q = q;
            KIDD = kidd;
            KGDD = kgdd;
            dt = dt_;
        }

        ~Source() // 析构函数
        {
        }

        // ========模型参数======== //
        public double SM;   // 表层自由水蓄水容量/mm ，敏感
        public double EX;   // 表层自由水蓄水容量方次，不敏感，1.0~1.5
        public double KG;   // 表层自由水蓄水库对地下水的日出流系数，敏感
        public double KI;   // 表层自由水蓄水库对壤中流的日出流系数，敏感

        // ========模型状态======== //
        public double R;    // 总径流量，mm
        public double RS;   // 地面径流，mm
        public double RI;   // 壤中流，mm
        public double RG;   // 地下径流，mm
        public double PE;   // 净雨量，mm
        public double FR;   // 本时段产流面积比例
        public double S0;   // 本时段初的自由水蓄量，mm
        public double S;    // 本时段的自由水蓄量，mm

        public double M;    // 一天划分的计算时段数
        public double KID;  // 表层自由水蓄水库对壤中流的计算时段出流系数，敏感
        public double KGD;  // 表层自由水蓄水库对地下水的计算时段出流系数，敏感
        public double SMM;  // 全流域单点最大的自由水蓄水容量，mm
        public double SMMF; // 产流面积上最大一点的自由水蓄水容量，mm
        public double SMF;  // 产流面积上的平均自由水蓄水容量深，mm
        public double AU;   // 相应平均蓄水深的最大蓄水深，S0值对应的纵坐标，mm
        public double RSD;  // 计算步长地面径流，mm
        public double RID;  // 计算步长壤中流，mm
        public double RGD;  // 计算步长地下径流，mm
        public double FR0;  // 上一时段产流面积比例
        public int N;       // N 为计算时段分段数，每一段为计算步长
        public double Q;    // Q 是每个计算步长内的净雨量，mm
        public double KIDD; // 表层自由水蓄水库对壤中流的计算步长出流系数，敏感
        public double KGDD; // 表层自由水蓄水库对地下水的计算步长出流系数，敏感
        public double dt;   // 模型计算时段长,h
    }

}
